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Statistical  Control  Paradigm  for  Aerospace  Structures 
Under  Impulsive  Disturbances 

Khanh  D.  Pham*  and  Lawrence  M.  Robertsom 

Air  Force  Research  Laboratory,  Kirtland  AFB,  New  Mexico,  87117,  U.S.A. 

In  this  paper  the  newly  developed  statistical  control  theory  is  revisited  to  autonomously 
control  the  satellite  attitude  as  well  as  to  provide  a  means  of  actively  attenuating  impulsive 
disturbances  caused  by  servicing  dock  and  space  debris.  Simulations  are  performed  using 
several  docking  and  collision  scenarios.  The  simulation  results  indicate  that  the  existing 
attitude  control  system  with  an  innovative  and  robust  statistical  controller  design  shows 
significant  promise  for  use  in  attitude  hold  mode  operation  despite  the  presence  of  impulsive 
disturbances. 

Nomenclature 

Ftc  Angular  momentum  of  the  satellite  about  its  center  of  mass  measured  in  inertial  frame  {N} 

Hw  Angular  momentum  of  the  wheel  cluster 
Jc  Satellite  inertia 
Jw  Reaction  wheel  inertia 

to S/N  Angular  velocity  vector  in  fixed-body  reference  frame  {S'} 

aiw  Angular  velocity  of  the  reaction  wheels 

L  Wheel  orientation  matrix 

Tthrt  Absolute  torque  due  to  the  thrusters 

tw  Absolute  torque  due  to  the  reaction  wheels 

igrav  Torque  due  to  the  earth  gravity  gradient 

Tmag  Torque  due  to  the  earth  magnetic  field 

Taero  Torque  due  to  atmospheric  drag 

rsrp  Torque  due  to  solar  radiation  pressure 

I.  Introduction 

Recent  work1-5  reported  at  several  American  Control  Conferences  showed  that  the  statistical  control 
design  has  been  performing  quite  competitively  with  other  modern  control  techniques.  In  the  statistical 
control  formulation,  a  state-feedback  controller  is  designed  to  minimize  the  objective  function  representing 
a  linear  combination  of  finite  cumulant  indices  of  a  finite  horizon  integral  quadratic  performance  measure 
associated  to  a  linear  stochastic  system.  A  dynamic  programming  approach  is  used  to  obtain  the  optimal 
control  solution.  This  control  algorithm  is  then  applied  to  attenuate  dynamical  effects  due  to  impulsive 
disturbances.  Motivations  for  this  arise  from  the  need  to  study  control  problems  related  to  antennas  on  the 
space  station  subject  to  impact  from  space  debris  and  active  damping  of  vibrations  of  flexible  structures 
caused  by  impact  forces.  The  outline  of  the  paper  goes  as  follows.  First,  a  brief  modeling  of  satellite  with 
reaction  wheels  and  thrusters  will  be  given  below.  Then  some  stabilizing  controllers  minimizing  the  first 
three  cost  cumulant  indices  from  space  disturbances  to  the  controlled  output  are  subsequently  designed. 
Finally,  the  results  of  simulation  are  presented  along  with  some  discussions. 

*  Aerospace  Engineer,  Space  Vehicles  Directorate,  3550  Aberdeen  Ave  SE,  and  AIAA  Member. 

^ Dynamics  &;  Controls  Program  Manager,  Space  Vehicles  Directorate,  3550  Aberdeen  Ave  SE,  and  AIAA  Senior  Member. 

Copyright  ©  2007  by  the  American  Institute  of  Aeronautics  and  Astronautics,  Inc.  The  U.S.  Government  has  a  royalty-free 
license  to  exercise  all  rights  under  the  copyright  claimed  herein  for  Governmental  purposes.  All  other  rights  are  reserved  by  the 
copyright  owner. 
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II.  Mathematical  Model  of  Satellite  with  Reaction  Wheels  and  Thrusters 

This  section  gives  the  reader  an  overview  of  a  simple  mathematical  model  for  center  of  mass  steering  and 
attitude  control  of  a  Low-Earth-Orbiting  (LEO)  satellite  whose  detailed  definition  is  given  in  Won.6  This 
model  assumes  that  the  largest  disturbance  is  the  gravity  gradient  torque  while  other  space  disturbances 
are  small  and  can  thus  be  modeled  as  a  stationary  Wiener  random  process.  Moreover  the  internal  torques 
generated  by  a  cluster  of  four  reaction  wheels  and  three  thrusters  are  for  use  in  controlling  the  attitude  of 
the  satellite.  Then  the  angular  momentum  of  the  satellite  with  three  thrusters  and  a  cluster  of  four  reaction 
wheels  becomes 

Hc  =  Jcas/N  +  LtJwuw  (1) 

where  Jw  is  the  mass  moments  of  inertia  of  the  reaction  wheels;  u)w  is  the  angular  velocity  of  the  wheels; 
and  L  is  the  wheel  orientation  matrix 

0  0  0  1  I  cos  a  sin  (3  sin  a  sin  / 3  cos  (3 

0  0  0  _  —  sin  a  sin  / 3  cos  a  sin  f3  cos  [3 

0  0  J£j  0  —  cos  a  sin  (3  —  sin  a  sin  (3  cos  (3 

0  0  0  J44  J  |_  sin  a  sin  (3  —  cos  a  sin  (3  cos/3 

Furthermore,  the  transpose  of  the  wheel  orientation  matrix  L  is  denoted  by 

'  Lt(  1,1)  Lt(  1,2)  Lt(  1,3)  Lt(l,4)  ' 

LT  =  Lt(  2,1)  Lt(  2,2)  Lt(  2,3)  Lt(  2,4) 

_  Lt( 3,1)  Lt( 3,2)  Lt{ 3,3)  Lt( 3,4)  _ 

with  a  =  45  degrees  and  (3  =  54.74  degrees.  By  taking  the  inertial  time  derivative  of  (1),  the  nonlinear 
satellite  attitude  dynamics  can  be  described  in  terms  of  vector  and  matrix  notations 

Jc3js/n  +  LT Jw3jw  +  us/N  (^Jclos/n  +  LT Jwun  =  rthrt  +  rgrav  +  w,  (3) 

where  the  stationary  Wiener  process  w  and  the  cross  product  matrix  Cbs^N  are  defined  as 

0  —  w3  w2 

A  1  1  ~S/N  A  „ 

W  —  T'aero  i  Tmag  +  Tsr p  ,  id  —  CJ3  (J  — UJi 

—10  2  iO\  0 

Furthermore,  the  angular  momentum  of  the  wheel  cluster  H\y  and  the  absolute  torque  due  to  the  reaction 


wheels  Tw  are  related  by 

Hw  =  tw  ,  (4) 

Hw  =  Jwcuw  +  JwLcus/n  .  (5) 

From  the  equation  (5),  it  is  possible  to  write 

vw  =  J^rw  -  Lus'N  ,  (6) 

which  results  in 

(Jc  -  LtJwL)  ios/N  =  -Cos'N  (jcus>N  +  LtJwuw)  -  Lttw  +  rthrt  +  rgrav  +  w.  (7) 

Let  Ja  —  Jc  ~  LT  JWL.  The  equations  (6)  and  (7)  can  then  be  rewritten  as  follows 
us/N  =  -Ja^S/N  (jCus,N  +  LtJwuw )  -  Ja1Lttw  +  J^Tthrt  +  J^Vav  +  J^w ,  (8) 


LOW  =  Jy^Tw  +  LJaOJS/N  ( JC0JS/N  +  LT  JWUJW^j  +  LJ LT  TW  -  I/J^Tthrt  ~  LJ^T grav  -  LJ^W  . 

(9) 
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Moreover,  the  orientation  of  the  body-fixed  reference  frame  S  relative  to  the  orbital  reference  frame  O  can 
also  be  described  by  introducing  the  time  dependence  of  quaternions  q  =  (91,  92,  93,  qi)T ■  Note  that  the 
orbital  reference  frame  O  =  {Ox,Oy,Oz}  is  rotating  about  the  Oy  axis  with  respect  to  the  Newtonian 
inertial  reference  frame  N  at  the  orbital  rate  to 0.  The  axes  of  O  are  chosen  such  that  the  roll  axis  Ox  is  in 
the  flight  direction,  the  pitch  axis  Oy  is  perpendicular  to  the  orbital  plane  in  the  negative  direction  and  the 
yaw  axis  Oz  points  from  the  satellite  to  the  earth  center.  The  benefits  of  using  quaternions  as  opposed  to 
Euler  angles  include  no  inherent  geometric  singularity  and  suitability  of  onboard  real-time  computation.  It 
is  shown  from  Wie10  that  the  kinematic  differential  equations  for  quaternions  can  be  shown  as  follows 


9i 

0  U)z  -UJy  LUX 

9i 

92 

1 

-cuz  0  LUX  UJy 

92 

93 

”  2 

LUy  -LUX  0  U)z 

93 

94  . 

-Wj  -UJy  - LOz  0 

.  94  . 

(10) 


where  ojs^°  =  (ijjx,LOy,LOz)T  is  the  rotation  rate  of  the  body-fixed  reference  frame  S  with  respect  to  the  orbit 
reference  frame  O,  expressed  in  the  body-fixed  reference  frame.  Moreover,  the  direction  cosine  matrix,  Cs^° 
is  needed  for  a  coordinate  transformation  that  maps  the  angular  velocity  from  u>s^°  to  lus^n 


Thus,  when 


ws/iv  =  ws/o  +  cS/O^O/N  ; 


c3/° 


Q1-Q2-Q3  +  94 

2(91 92  ^  9394) 
2(9i93  +  9294) 


2(9i92  +  9394) 
92  -  9i  -  93  +  94 
2(9293  ^  9i94) 


2(9i93  ^  9294) 
2(9293  +  9i94) 
93  -  Qi  ~  92  +  94 


w°'N  = 


0 

-w0 

0 


(11) 

(12) 


(13) 


is  the  inertial  rotation  rate  of  the  orbital  reference  frame  with  respect  to  the  inertial  frame,  the  uix,  u>y,  and 
ojz  in  (10)  can  be  replaced  according  to  the  following  equations 


Ux  =  wi  +  2w0(9i92  +  9394) ,  (14) 

vy  =  w2  +  w0(92  -  9i  -  93  +  94) .  (15) 

wz  =  w3  +  2w0(9293  -  9i94)  ■  (16) 

The  gravity  gradient  torque  about  the  satellite’s  mass  center  is  given  by  Wie10 

^gravity  —  3w>q Oz  X  JcOz  ,  (If) 


where  the  vector  Oz,  represents  the  projection  of  the  orbit  radius  vector  onto  the  body-fixed  reference  frame, 

S 


Thus,  the  equation  (17)  becomes 


oz  =  cU°h  +  cU°h  +  cl 'l°h  ■ 


(18) 


0 

nsio 

°33 

ns/o  j 

L/23 

^gravity  —  3^0 

rs/o 

L'33 

0 

ns/o 

L  13 

ns/o 

'-'23 

rs/o 

L  13 

0 

0 

—  3  uJq 

(93  -  1l 

-  92  +  9l) 

—  2(9293  +  9194) 

Jll 

0 

0 

r  nsl°  1 

c'13 

0 

+22 

0 

rs/o 

c'23 

0 

0 

+33 

rs/o 

.  L'33 

-(93-9i  -92+91)  2(9293  +  9i94) 

0  -2(91 93  -  9294) 


2(9i93  -  9294)  0 


+11 

0 

0 

2(9193  -  9294) 

0 

+22 

0 

2(9293  +  9194) 

0 

0 

+33 

.  (93  -  9i  -  92  +  9l)  . 

(19) 
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Finally  the  the  rotational  motions  of  the  satellite  with  three  rotational  degrees  of  freedom  can  be  described 
by  the  kinematic  differential  equations 


91  =  -  (W154  -  w2g3  +  u3q2  +  w 0g3)  ,  (20) 

92  =  ^  (wig3  +  w294  -  u3q  1  +  W094)  ,  (21) 

93  =  7^  (-W192  +  w29i  +  w394  -  w09i)  ,  (22) 

94  =  -  (— <x>i9i  —  w292  —  w393  —  W092)  i  (23) 

JauS/N  =  - uS/N  (jcuS/N  +  LTJwcow)  -  Lttw  +  rth rt  +  xgrav  +  w  ,  (24) 


W^  =  LJ^CoS/N  (jCUS,N  +  LTJWUW )  +  +  J^1)  TW  -  LJ^Tthrt  -  LJ^Vgrav  -  . 

(25) 

which  can  be  further  rewritten  compactly  using  the  mapping  /  :  R11  xl'  x  R  1— >  R11  whose  rule  of  action 
is  given  by 

x  =  f(x,u,t)  +  v ,  te[t0,tf]  (26) 

wherein  the  state  variable,  x  £  R11;  control  variable  constrained  to  some  compact  set,  u  £  R7;  and  combined 
torque  disturbances,  v  £  R3  are  defined  by 


Xl 

9i 

*2 

92 

*3 

93 

Ul 

Tw(  1) 

X4 

94 

U2 

rw(  2) 

x5 

Wi 

U3 

Tw(  3) 

C*4x3 

x6 

A_ 

w2 

,  u  = 

U4 

A_ 

Tw(  4) 

,  v  =  Gw  = 

Ja 

x7 

w3 

u5 

Tthrt(l) 

if 

_ 1 

*8 

,  ,w 

W 1 

■u6 

!"thrt  (2) 

Xg 

,  ,w 
w2 

.  u7  _ 

Xthrt  (3) 

X10 

,  ,w 
w3 

*11  - 

,  ,w 

L  w>4  J 

and  to  and  tf  are  the  beginning  and  ending  mission  lifetimes.  Since  the  high  fidelity  nonlinear  model  in 
(20)-(25)  is  too  complicated  for  use  in  control  design.  It  is  necessary  to  develop  a  simpler  linearized  model 
which  can  approximate  the  system  over  a  wide  range  of  conditions.  The  following  linearized  model  assumes 
a  circular  orbit,  and  linearization  is  performed  about  the  equilibrium  nadir-pointing  attitude.  Clearly,  it  is 
easy  to  see  that  there  exist  two  trivial  equilibrium  points  ( x*,u *)  for  the  deterministic  system  x  =  f(x,u,t) 


x *  =  (0,0, 0,1,0, -w0, 0,0, 0,0,0), 

u*  =  (0,0, 0,0, 0,0,0) , 


x*  =  (0,0,0, -1,0, -wo,  0,0, 0,0,0), 
u*  =  (0,0, 0,0, 0,0,0) , 


such  that 


f(x*,u*,t)=  0,  t£[t0,tf]. 


Once  the  satellite  is  in  orbit,  the  satellite  is  controlled  to  keep  its  pointing  direction  about  a  pre-computed 
orientation.  In  this  case,  a  reference  orientation  is  defined  as  a  state  x*  in  which  the  deviation  of  orientation 
is  preferably  zero  under  a  steady  state  condition.  The  control  objective  is  to  regulate  the  orientation  of 
the  satellite  about  these  desired  orientations.  The  satellite’s  tracking  of  itself  to  the  prescribed  orientation 
is  equivalent  to  the  expression  x  —  x*  being  zero.  This  indicates  that  linearization  of  (26)  about  some 
reference  orientation  can  be  useful  for  studying  the  effect  of  perturbation  away  from  the  desired  orientation 
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and  remains  accurate  for  sizable  changes  in  the  satellite’s  angular  velocity.  Reaction  wheels  and  electric 
thrusters  can  be  used  to  actively  drive  observed  perturbations  back  to  zero.  The  expansion  of  the  right-hand 
side  of  the  equation  (26)  in  a  first-order  Taylor  series  expansion  in  terms  of  small  perturbations  proceeds 
with 


Sx  =  x  —  x*  , 
Su  =  u  —  u*  , 


where  the  reference  orientation  and  the  nominal  control  input  are 


sk 

X2 

sk 

^3 

< 

x% 

U*2 

sk 

£5 

sk 

u3 

sk 

^6 

* 

,  u  = 

u\ 

Xj 

sk 

u5 

00  * 

K 

sk 

^9 

sk 

w7 

* 

XW 

sk 

^11  J 

The  linearization  of  (20)-(25)  about  the  reference  orientation  yields 

Sx(t)  =  A(x*(t),u*(t))6x(t)  +  B(x*(t),u*(t))Su(t)  +  v(t) ,  Sx(to)  =  xo  ,  (27) 

where  matrix  coefficients  A(-,  •),  and  B(-,-)  are  given  by 


where 


A, 


A(x*(t),u*(t)) 


df 


dx 


(x*  ,U*) 


B(x*(t),u*(t)) 


df 


du 


(x*  ,U*) 


l:4*£l:4 

^1:40:5:7 

^1:4^8:11 

-4*5 :  7*1 :4 

-4*5: 7*5 : 7 

*4*5;  7*8:  H 

*4*8:11*1:4 

-4*8:11*5:7 

-4*8:  H  *8:  H 

-®*1:4«1:4 

*3*1:4115:7 

-®*5:7«1:4 

***5:7115:7 

-®*8:11«1:4 

4**8:11115:7 

(28) 


(29) 


0 

0 

w0 

0  " 

'  1 

2 

0 

0  " 

'  0 

0 

0 

0  " 

0 

0 

0 

0 

4 

5  -^Xi:4X5:7  ~ 

0 

1 

0 

4 

5  ^3:1:4^8:11  _ 
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0 

0 

0 

—w0 

0 

0 

0 

0 

2 

0 

1 

2 

0 

0 

0 

0 

0 

0 

0 

0  _ 

_  0 

0 

0  _ 

_  0 

0 

0 

0  _ 

A, 


( Jc  ~  LT  JWL)  1 


6wg(J33  -  J22)  0  0  0 

0  6wg(J33  -  Jn)  0  0 

0  0  0  0 


Ax 5;7X5;7  —  (Jc  J  Jw  J) 


0  0  —  Wo(<^22  —  J33) 

0  0  0 

-ixo(Jn  —  J22)  0  0 


-^•0:5:7X8:11  (Jc  ~  LT JWL)  1 


w0Lt(3,l)Jff 

0 


uoLt(3,2)jW  w0Lt(3, 3)  J33 
0  0 
-u0Lt(l,2)J%  -w0it(l,3)J^ 


w0Lt(3,4)jjX 

0 

~ 1 ^oi't(l)  4)  J44 
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^4x8:11X1:4  ^ 1  ^4x8:11X5:7  ^Xs, 7X5 ,7  I  ^4x8:11X8:11  ^ ^-Xf, ; 7 X8: 1 1  > 


Ij, 


'  0 

0 

0 

0  ' 

'  0 

0 

0  ' 

0 

0 

0 

0 

)  -®Xl:4«5:7  - 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

_  0 

0 

0 

0  _ 

_  0 

0 

0  _ 

BX5;Tu1:4  =  -( Jc  -  LT  JWL )~1LT  . 


B 


-Sx8:llHl:4  =  L{Jq  ~  L  J W  L)  L  +  J, 


/-1 
'iV  5 


X5:7U5:7  —  {Jc  —  LT  JwL)  1 , 

1x5 ; 7  =  —L(Jc  -  LT JwL)-1 . 


B 


III.  Model  of  Satellite  In  Events  of  Servicing  Dock  and  Debris  Impacts 

Next,  the  conservation  of  angular  momentum  is  utilized  in  the  event  of  docking  and  space  debris  impacts. 
The  assumptions  made  herein  include  not  limited  to  the  case  where  a  servicing  satellite  that  is  slowly 
attaching  to  the  satellite,  or  the  case  where  high  speed  debris  is  embedded  in  a  perfectly  inelastic  collision. 
Thus,  the  angular  momentum  immediately  before  and  after  the  impact  event  must  match 


Hsat  T  LI Docking / DebriS 


L^Sat  T  H Docking  /  Debris 


J  Post— Impact 


(30) 


J  Pre— Impact 

It  is  well  to  note  that  the  angular  momentum  of  the  satellite  before  and  after  the  impact  are 

Lt Sat\pre— Impact  0? 

Hsat  |  Post-Impact  =  J Sat+Docking  /  Debris^  ^  1 

and  the  angular  momentum  of  either  a  servicing  satellite  or  space  debris  is  given  by 

H Docking  /  Debris\  Pre- Impact  T  Docking  /  Debris  ^  ^klDocking  /  Debris'^  Docking  /  Debris  1 

H Docking  /  Debrisl  Post- Impact  0  ; 

wherein  r Docking / Debris  is  the  position  vector  from  the  satellite’s  center  of  mass  to  the  point  of  impact 
and  mDocking / Debris^ Docking / Debris  is  the  linear  momentum  of  either  a  servicing  satellite  or  space  debris. 
Therefore,  the  equality  (30)  can  be  readily  reduced  to 


I" Docking / Debris  ^  Docking/ Debris^  Docking  /  Debris  J Sat.+ Docking / Debris^I  ^  : 


(31) 


provided  Jsat+Docking/ Debris  =  Jc  +  A  Jc  is  the  combined  moment  of  inertia  of  the  satellite,  servicing 
satellite  or  space  debris.  Resolving  f Docking /Debris  and  v Dockingi /Debris  hrto  the  body-fixed  reference  frame  S 

C Docking / Debris  =  ^*1^1  T  T2S2  T  I"3^3  j 
V Docking / Debris  —  ^1^1  T  V2&2  T  ^3^3  5 

yields  the  angular  velocities  of  the  satellite  in  the  event  of  either  docking  or  collision 

r2V3  -  r3v2 


m  Docking / Debris  jll 

°  Sat-\- Docking  /  Debris 

_  rxv 3  -  r3^i 

Cc?2  in  Docking / Debris  t22 

0  Sat-\-  Docking  /  Debris 

r\V2  —  r2v  1 

CJ3  m  Docking j  Debris  t33 

°  Sat-\- Docking  /  Debris 


(32) 

(33) 

(34) 


From  the  equation  (31),  it  is  observed  that  a  soft  docking  would  be  possible  if  the  servicing  satellite  had 
small  mass,  docked  slowly  and  was  close  to  the  center  of  mass  of  the  satellite.  For  simplicity,  impact  forces 
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caused  by  either  the  docking  event  or  space  debris  are  modeled  as  finite  duration  events  and  it  is  further 
fruitful  to  assume  that  each  impact  duration  is  very  short  and  is  denoted  as  ti+ 1  —  tl  =  hi.  Thus,  the  impact 
can  be  treated  as  an  impulsive  event  and  the  impact  force  wd(ti)S(t  —  tt)  generated  by  the  ith  impact  and  the 
resulting  angular  velocity  of  the  satellite  (U),  u>2 (tj),  and  u>3 (ti)  caused  by  the  ith  impact  force  is  defined 
as 


m  Docking  /  Debris(ti) 


r2(ti)v3(ti)  -  r3(ti)v2(ti) 
I11  (t  )  ’ 

Sat-\-  Servicing  /  Debris''  1 ' 


^2  yj'i)  —  "‘'Docking/ Debris  {"i)  j22 


1  Sat-\-  Servicing  /  Debris 


( U )  ’ 


,  ,  ,  ,ri/ti)v2{ti)  -  r2(ti)v1(ti) 

^3  V'i)  "’'Docking/ Debris  V'i) 


'  S  at-\-  S  ervicing  /  Debris 


( U )  ’ 


and  the  resultant  change  of  the  quaternions  qi(ti),  q2(U),  q3(ti),  and  q±(ti)  are  obtained  via  the  kinematic 
equations  (10).  It  is  worth  noting  that  the  effect  of  impulsive  disturbances  causes  possible  jumps  of  state 
variables.  The  state  space  realization  of  the  satellite  in  the  event  of  the  docking  or  space  debris  impacts  is 
therefore  a  hybrid  system  which  contains  both  continuous-time  and  discrete-time  components  and  can  be 
expressed  as  follows, 


Sx(t)  =  A(t)8x(t)  +  B(t)5u(t)  +  G(t)w(t)  +  Ed(ti)wd(ti)8(t  —  ti) ,  8x(to)  =  Xq  ,  (35) 

*3=1 

where  8{t)  is  the  Dirac  Delta  function;  the  impulsive  disturbances  wd(U)  €  h([to,tf]  ,  R3)  the  Hilbert  space 
of  square-summable  Revalued  sequences  with  the  norm  defined  by  jliCdH3,  =  X)t  e[t0  tf]  wd(li)wd(U)i  the 
states  Sx(t)  G  C([to,t f}'^11))  the  subset  of  Hilbert  space  of  R11-valued  square-integrable  process 

on  [to,tf]  that  are  adapted  to  the  cr-field  Tt  generated  by  w(t )  with  the  norm  defined  by  H&rHf,  = 
E  {fto  $xT  (t)Sx(t)cIt^  <  oo;  and  the  control  input  8u(t)  G  L^t(H;C([to,t/];R7))  the  subset  of  Hilbert 
space  of  Revalued  square  integrable  process  on  [t0 ,tf]  that  are  adapted  to  the  u-field  Tt  generated  by  w(t) 
with  the  norm  defined  by  ||<5w|||2  =  E  j  f/f  SuT (r)<5M(r)(ir|  <  oo.  The  impact  instants  ti  G  [t0,tf]-  No¬ 
tice  that  continuous-time  bounded  coefficients  A(t),B(t),  and  G(t)  are  defined  as  before  with  appropriate 
dimensions  whereas  the  discrete-time  matrix  coefficient  Ed(ti)  is  given  by 


Ed(U) 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
10  0 
0  10 
0  0  1 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


(36) 


It  is  also  interesting  to  represent  the  hybrid  system  (35)  in  the  form  of 

8x(t)  =  A(t)Sx(t)  +  B(t)5u(t)  +  G(t)w(t) ,  t  yf  t j ,  8x(to)  =  x3  ,  (37) 

Sx(t+)  =  8x(t)  +  Ed{ti)wd{ti) ,  t  =  U  ,  *  =  1,2,...,  (38) 

where  Sx(t/ )  are  the  values  of  state  dynamics  immediately  after  the  impulses,  and  the  time  limit  t+  is  defined 
as  lime^0(i*  +  e)  for  e  G  R+.  It  is  important  to  note  that  the  state  variables  Sx(t)  are  right  continuous  and 


7  of  17 


American  Institute  of  Aeronautics  and  Astronautics  Paper  AIAA-2007-1755 


may  be  left  discontinuous  due  to  the  discrete  time  jumps.  The  continuous  and  discrete  system  parts  described 
by  (37)- (38)  are  rewritten  as 

Gc  '■  |  Sx(t)  =  A(t)5x(t)  +  B(t)Su(t)  +  G(t)w(t) ,  t  ,  Sx(t0)  =  xq  ,  (39) 

Gd  ■  |  Sx(t+)  =  Sx(t )  +  Ed(ti)wd(ti) ,  i  =  1,2,...,  (40) 

As  depicted  in  Figure  1,  a  reference  orientation  is  fed  into  the  satellite  controller  subsystem  which  then 
calculates  the  wheel  and  thruster  commands  based  on  the  target  orientation  and  the  current  orientation. 
Consequently,  it  results  in  a  change  in  reaction  wheel  angular  accelerations  and  thruster  torques  and  thus 
the  satellite  motion  is  altered.  The  motion  is  sensed  in  the  body  frame  of  reference  and  is  changed  into  an 
inertial  frame  by  the  kinematic  equations. 


Kinematic 

Gyro 

Inertial  Orientation 

Equations 

Body  Rates 

Figure  1.  Satellite  Functional  Block  Diagram 


IV.  Statistical  Control  Development 

As  the  satellite  is  launched  into  an  intended  orbit,  it  never  remains  in  the  ideal  orientation.  In  addition 
to  servicing  dock  and  space  debris  impacts,  the  external  forces  present  in  space  created  by  32,  solar  radiation 
pressure,  and  aerodynamic  torques  cause  perturbations  to  this  ideal  orientation.  Thus,  the  linear  time- 
varying  system  with  finite  discrete  jumps  can  now  be  rewritten 

Sc  '■  |  dx(t)  =  (A(t)x(t)  +  B(t)u(t))  dt  +  G(t)dw(t) ,  t^U,  x(to)  =  xo  ,  (41) 

Gd  ■  {  xd(t+ )  =  xd(t)  +  Ed(ti)wd(ti) ,  t  =  ti,  i  =  1,2,...,  (42) 

where  A  £  C([f0,  t /] ;  Rnxn),  B  £  C([t0,  t /] ;  Rnxm),  G  €  C([to,tf];  Rnxp)  and  the  system  noise  w(t)  £  Rp 
representing  all  external  forces  in  space  other  than  space  debris  collisions  and  servicing  docking,  is  a  stationary 
Wiener  process  with  correlation  of  increments  E  {[w(i)  —  tu(cr)]  [w(f)  —  w(ct)]t}  =  W\t  —  oj.  Also,  assume 
the  initial  state  x(to)  =  Xq  is  known.  For  a  given  (to,Xo),  associate  with  the  nominal  system  (41)  a  finite 
time  integral  quadratic  form  (IQF)  random  cost  J  :  C{\t0,tf\,  Rm )  i— >  R+  such  that 

rtf 

J{u)  =  xT(tf)Qfx(tf)  +  /  [xt(t)Q(t)x(t)  +  ut(t)R(t)u(t)]  dr  ,  (43) 

Jt0 

where  Qj  £  Rnxn,  Q  £  C([io,  £/];  Rnxn)  and  R  £  C([t0,tf};  Rmxm)  are  symmetric  and  positive  semidefinite 
with  R(t)  invertible. 

As  observed  by  Liberty,7  all  cumulants  of  any  IQF  cost  in  the  system  state  of  linear  dynamical  systems 
are  quadratic  affine  in  the  initial  state.  In  the  case  of  state  measurement,  it  is  thus  fruitful  to  assume  the 
control  input  u{t)  to  be  a  linear  time-varying  feedback  law  given  by 

u(t)  =  K(t)x(t) ,  te[t0,tf]  (44) 

in  which  K  £  C([to,  t/];  Rmxn)  is  an  admissible  feedback  gain  whose  definition  will  be  clear  shortly. 

For  a  given  initial  condition  (tg,Xo)  and  an  admissible  gain  I\,  the  kth  cost  cumulant  of  random  cost  J 
with  fixed  k  £  Z+  is  given  by  Liberty7  by 


Kk{t0,  xo',  K)  =  XoH{t0,  k)x o  +  D (t0,  k), 


(45) 
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where  cumulant-building  variables  H(a,k)  and  D{a,k )  evaluated  at  a  =  to  are  satisfying  the  differential 
equations  (suppressing  time  argument  of  matrix  coefficients) 


4-H(a,  1)  =  -(A  +  BI<)TH(a,  1)  -  H{a,  1){A  +  BK)  -  I\T RK  -  Q  ,  (46) 

da 

7  i  1  c\  •  | 

—H(a,i)  =  ~(A  +  BK)TH(a,  i)  -  H(a,i)(A  +  BK)  -V—  l'  ^H(a,j)GWGT H(a,i  -  j),  (47) 

da  “( j'(*  ~  J)- 

4-D(ot,i)  =  -Tr\H(a,i)GWGT\  ,  1  <i<k,  (48) 

da 

with  the  terminal  conditions  H(tf ,  1)  =  Qf,  H(tf,  i)  =  0  for  2  <  i  <  k,  and  D(tf,  i)  =  0  for  1  <  i  <  k. 

Now  it  is  convenient  to  denote  the  right  members  of  the  equations  (46)-(48)  by  the  symmetric  mappings 

Ti  (a,  H,  K)  =  -(A  +  BK)TH(a ,  1)  -  H(a,  1)(A  +  BK)  -  KT RK  -  Q , 

Fi(a,H,K)  =  -(A  +  BK)TH(a,i)  -  H(a,i)(A  +  BK)  Y  „H(a,j)GWGT  H(a,i  -  j), 

Gi(a,H)  =  -Tr{H(a,i)GWGT}  ,  1  <i<k, 

with  fc-tuple  matrix  variables 


H{a)  =  (H(a,l),...,H(a,k)), 

V(a)  =  (D(a,l),...,D(a,k)). 

Thus,  the  product  system  of  (46)-(48)  is  written  as  follows 

^W(a)  =  F(a,H,K),  H(tf)=Hf, 

A-V(a)  =  Q(a,H),  V(tf)=Vf, 

where  hi  /  =  (Qf,  0, . . . ,  0)  and  Df  =  (0, . . . ,  0).  Next  denote  explicitly  the  dependence  of  H(a)  and  V(a)  on 
the  admissible  control  gain  K  by  H(a,K)  and  V(a,K).  Hence,  the  performance  index  in  the  kCC  control 
problem  follows. 

Definition  1  (Performance  Index) 

Fix  a  k  £  Z+  and  a  sequence  p  =  {pi  >  0}^=1  with  p i  >  0.  With  (to,Xo)  given,  the  performance  index  of 
the  state  feedback  kCC  control  problem  is  defined  by 

k 

4>0  (t0,H(t0,K),V(t0,K))  =  y^jiiKi(to,x0;K) 

2=1 

k  k 

=  xo^2  K)x0  +  ^2  K) ,  (49) 

2—1  2=1 

where  the  real  constants  pi  give  parametric  control  freedom. 

The  class  K.^  of  admissible  feedback  gains  is  then  defined. 

Definition  2  (Admissible  Feedback  Gains ) 

Let  the  compact  subset  K  C  Rmxn  be  the  allowable  set  of  gain  values.  For  a  given  k  £  Z+  and  a  sequence 
k-  =  {ft  >  0}*_i  with  pi  >  0,  let  K.^  be  the  class  of  C([to,tf}\ Rmxn)  with  values  K(-)  £  K  for  which  there 
exist  solutions  to  the  dynamic  equations  of  motion 


—'H(a)  =  T(a,  7i(a),K(a)) ,  B(tf)  =  Hf  , 

(50) 

-^D(a)=Q(a,U(a)),  V(tf)=Vf. 

(51) 
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The  fcCC  control  optimization  problem  is  now  stated. 

Definition  3  (Optimization  Problem) 

Suppose  both  k  £  Z+  and  the  sequence  p  =  {pi  >  0}^=1  with  pi  >  0  are  fixed.  Then  the  finite  horizon  state 
feedback  kCC  control  optimization  may  be  defined  by 


min  (j)0  (t0,  H(t0,  K),T>(to,  K))  ,  (52) 

G/v S  [A 

subject  to  the  dynamic  equations  of  motion,  for  all  a  £  T , 

^H(a)  =  f(a,H(a),K(a)),  H(tf)=Hf,  (53) 

da 

^V(a)  =  g(a,H(a)) ,  V{tf)  =  Vf.  (54) 

This  “Mayer  form”  optimization  problem  can  be  solved  by  applying  an  adaptation  of  the  HJB  equation  and 
verification  theorem  of  dynamic  programming  given  in  Fleming  and  Rishel.8  Readers  who  are  interested  in 
seeing  the  detailed  adaptation  and  associated  proofs  should  refer  to  the  reference.9 

Theorem  1  (HJB  Equation) 

If  there  exists  an  optimal  control  gain  K*  £  K fi,  then  the  partial  differential  equation  of  dynamic  programming 


min_ 

KGI< 


!v(^,Z) 


d 

dvec(y ) 


V{e,y,Z)vec(P(e,y,K)) 


d 


dvec(Z) 


V(e,y,Z)vec(g(e,y))  >  =o. 


(55) 


is  satisfied.  Here  V(-,-,  •)  is  the  value  function. 

The  verification  theorem  stated  in  the  notation  used  here  is 


Theorem  2  (Verification  Theorem) 

Fix  a  k  £  Z+  and  the  sequence  p  =  {pi  >  0}f=1  with  p\  >  0.  Also,  let  the  candidate  value  function 
W(e,y,Z)  be  a  continuously  differentiable  solution  of  the  HJB  equation  (55)  which  satisfies  the  boundary 
condition 

W(t0,H0,V0)  =  <t>o(t0,Ho,V0)  .  (56) 

Let  K  £  K. ^  be  any  admissible  control  gain  and  let  H(-)  and  V(-)  be  the  corresponding  solutions  of  the 
equations  of  motion  (53)-(5j).  Then  W(a1‘H(a),V(a))  is  a  non-increasing  function  of  a.  If  K*  £  is  a 
control  gain  with  corresponding  solution,  TL*(-)  and  V*(-)  of  the  equations  (53)-(5f)  such  that  for  a  £  [to,  f/] 


j-W(a,  H*(a),V*(a))  +  .  W(a,  H*(a),V*(a))vec(T(a,  H* (a) ,  K*  (a))) 
oe  ovec(y) 

+  dJc{z )  W(a,  H*  (a) ,  V*  ( a))vec(g(a ,  TL*  (a)))  =  0,  (57) 

then  K*  is  an  optimal  control  gain  in  and 

W(e,y,Z)  =  V(e,y,Z).  (58) 

where  V(e,y,  Z)  is  the  value  function. 

For  any  e  £  [t0,t/],  the  cumulant-generating  states  of  the  product  system  (53)- (54)  defined  on  the  interval 
[to,e]  have  terminal  values  denoted  by  U(e)  =  y  and  11(e)  =  Z.  Observe  that  the  performance  index  (49)  is 
quadratic  affine  in  terms  of  arbitrarily  fixed  Xq,  a  hypothesis  solution  to  the  HJB  equation  (55)  has  the  form 

k  k 

W(e,  y,  Z)  =  xlYJ  +  ^(£))*o  +  y  +  %(e)) . 

i=l  i= 1 


where  £i  G  C([to,tf];  Rnxn)  and  %  G  C([to,tf];  R)  are  unknown  functions  of  time.  Placing  this  guess  into 
(55)  yields 


miry 

KeK 


k  d  k  d  k  & 

y  p,  -j^£i(e)x o  +  y  Pi—Tfis)  +  Xq  y  PiFfie,  y,  K)x0  +  y  PiQi(e,  (F) 


i=  1 


i= 1 


'  de 


i= 1 


i=l 


=  0. 


(59) 
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Differentiating  the  bracket  expression  in  (59)  with  respect  to  K  yields  the  necessary  condition  for  an  ex¬ 
tremum  of  the  performance  index  (49)  on  [to,tf], 


K(e,y)  =  -R-1{E)BT{e)YJVryr, 


(60) 


r— 1 


where  Jir  =  ^i/fi i  and  ji-\  >  0.  Applying  the  feedback  gain  (60)  along  the  solution  trajectories  of  (53)-(54) 
yields 

d  k 

=  -Ar{e)Hi(e)  -  TCi(e)A(e)  -  Q(s )  +  Hi{e)B(e)R~~1(e)Bt (g)  Y  £s7i:s(g) 


S=  1 


+  ^ /xr7Yr(e)i3(£)i?  1(e)Bt(e)Hi(£)  —  '^2Tlr'Hr(£)B(e)R  1(e)Bt(e)  y^jisHs(£) ,  (61) 


r=l 


r— 1 


S  —  1 


=  -AT(£)Hi(e)  -  Hi(£)A(s)  +  TLi(£)B(e)R  l(e)BT(e)  ^  j2sHs(e) 


i- 1 


2i! 


+  V'/xr7fr(e)B(£)i?  1(e)BT(£)Wi(e)  -  V]  1  7ij(g)G(g)VrGT(g)7^-j(e) , 


cfe 


Z>i(e)  =  -Tr^g^g^G^g)}  . 


(62) 

(63) 


Replacing  (60)  into  the  bracket  of  (59)  yields  the  minimum 


Mi  Bi(s)  +  Me)A(e)  +  MiQ(^) 


,i=l 


i=l 


k 


-  ^)ir7ir(g)H(g)i?  1(£)BT(e)'^jni'Hi{£)  -  ’s^iHi'Hi{E)B{E)R  1(e)BT{£)  Y, %Hs(e) 


r= 1 


i=l 


i=l 
k  2—1 


AC  AC  AC  2—1  2  -| 

+  /X1  ^2’p,r'Hr{s)B(£)R^1  (e)BT  (£)^2fisHs{£)  +  r(i i 

r=l  s=l  i=2  .7  =  1  ^ ' 


s=l 
T I 


1  =  1 

ftTr{Wj(e)G(e)^Gr(£)}  .  (64) 

i=l  i=l 

Having  examined  (64),  £)(•)  and  7j(-)  are  chosen  to  satisfy  the  differential  equations 

d  k 

~r £i(e)  =  AT (£)Hi(e)  +  H\{e)A(e)  +  Q(e)  -  Hi(e)B(e)R~1(e)Bt(e)  Y  ^sHs{e) 


U- j)! 

k 


~Me)G{e)WG1  WHi-jie) 


x0 


d£ 


8=  1 


-  '^2j2rHr{E)B(£)R.  x{e)Bt  {e)Hi{e)  +  Y^rHr{£)B{£)R  1(g)HT(g)  ^  ^1sHs{e) 


(65) 


r=l 


r=l 


5  =  1 


— fi(s)  =  i4T(€r)Wi(e)  +  Hi(e)A(e)  -  Hi[e)B[e)R  1(e)BT{e)  ^fisHaje) 


2-1 


2z! 


~Y,V’r'Hr(e)B(£)R  1(£)Bt{£)U,{£)  +  Y  — - — —  Wj  (g)  G(g)fHGT(g)7f »_  y  (g) 

=1  Jj- 


> 


d£ 


%{e)  =  Tr  |71i(g)G(g)H/GT(g)}  ,  1  <  *  <  fc. 


(66) 

(67) 


At  the  boundary  condition  (56),  one  has 

k  k  k  k 

x0  y  ^  MiC^iO  “1“  ^i(^o))  ^0  H"  y  ^  Mi iO  0))  —  Xq  y  ^  Mi^iO  ^0  “1“  y  ^  Mi-^iO  • 


i=l 


i=l 


i=l 


i=l 
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Thus,  the  initial  conditions  are  found  to  be  £t(to)  =  0  and  7)  (to)  =  0  for  1  <  i  <  k.  Enforcing  the  initial 
conditions  and  comparing  (61)-(63)  to  (65)-(67)  uniquely  specify 

£»(e)  =  Hi{to)  -  Ki{e) ,  %(e)  =  Vffto)  -  Vfie) , 

for  any  £  £  [to,tf]  and  yields  the  value  function 

k  k 

yV(fo,W(to),^(io))  =  Xq  ^2  iMHi(t0)x0  +  y^/UjPj^o) , 

2=1  2=1 

which  satisfies  the  sufficient  condition  (57)  of  Theorem  2. 

Theorem  3  (Finite  Horizon  kCC  Control  Solution) 

Suppose  k  £  Z+  and  the  sequence  p  =  {pi  >  0}f=1  with  pi  >  0  are  fixed.  Then  the  optimal  state  feedback 
gain  K*  that  minimizes  the  performance  index  4> o  (to,  7t(to,  K),  V(to,  K))  is  given  by 

k 

K*(a)  =  -R-1(a)BT(a)^prH*(a,r),  a£[t0,tf]  (68) 

r=  1 

where  the  scalar,  real  constants  pr  =  Pi/p i  represent  parametric  degrees  of  design  freedom  and  solutions 
H(-,r)  satisfy  the  backward-in- time  matrix  differential  equations 

1)  =  -[A(a)  +  B(a)K*(a)]T H*(a,  1)  -  H*(a,  l)[A(a)  +  B{a)K*{a)\ 

-Q{a)-K*T(a)R(a)K*{a),  (69) 

—  H*(a,r)  =  —  [A(a)  +  B(a)K*(a)]TH*(a,r)  —  H*  (a,r)[A(a)  +  B(a)K*  (a)] 

r~1  q  r' 

V  — — -H*(a,s)G(a)WGT(a)H*(a,r  ~  s) ,  2  <r<k  (70) 

'  s\(r  —  s)! 

S—l  x  ' 

with  H*(tf,  1)  =  Qf  and  H*(tf,r)  =  0  for  2  <r  <k. 

V.  Simulation  Results 

This  section  contains  several  numerical  simulations  for  the  statistical  controllers  described  in  the  previous 
section.  The  low  earth  orbiting  satellite  is  sun-synchronous  remote  sensing  satellite  with  an  inclination  of 
98.13  degrees,  an  altitude  of  685  km  and  a  total  weight  of  509  kg.  It  is  assumed  that  the  angular  rate  and 
the  reaction  wheel  speeds  are  available  from  the  rate  gyro  and  the  reaction  wheel  tachometers,  respectively. 
The  orbital  rate,  u>o  is  0.0010636  rad/s,  the  total  moment  of  inertia  for  the  spacecraft  body  is 

’  294.62  0  0 

Jc  =  0  129.56  0  kgm 2 . 

0  0  209.76 

This  data  shows  that  the  x-axis  moment  of  inertia  is  the  largest,  followed  by  that  about  the  z-axis.  The 
moment  of  inertia  about  the  y-axis  is  the  smallest.  Hence  for  gravity  gradient  stability,  the  principal  x-axis 
of  the  body  should  be  aligned  with  the  orbit  normal  (pitch),  the  z-axis  should  be  aligned  with  the  velocity 
vector,  and  the  y-axis  should  be  the  yaw  axis.  The  body  axes  and  the  principal  axes  are  assumed  to  align 
with  each  other  as  the  products  of  inertia  are  negligible.  The  moment  of  inertia  of  the  reaction  wheels  is 

’  1  0  0  0  " 

Jw  =  0.01044  0  1  0  0  kgm2  , 

0  0  10 

_  0  0  0  1  _ 
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and  the  parameters  for  the  wheel  orientation  matrix,  L  are  given  as  a  =  45  degrees,  (3  =  54.74  degrees. 
Furthermore,  the  controller  design  parameters  are  selected  as  follows 

Q  =  10-/nxn,  i?=100-/7x7,  W  =  0.0045  •  /3x3  • 

Having  been  introduced  in  Section  III,  the  statistical  control  design  equations  (69)-(70),  of  course,  are  nonlin¬ 
ear,  and  so  represent  an  interesting  study  in  themselves.  In  order  to  assess  the  trends  of  performance  offered 
by  the  use  of  second  and  third  cost  cumulants,  it  was  chosen  in  the  present  instance  as  an  approximation 
approach  to  the  application  of  the  equations  (69)-(70).  The  time-invariant  statistical  controller  has  the  form 

I<  =  —R~1Bt  [Hx  +  ll2H2  +  ^3H3]  , 

where  H i,  H2,  and  H3  are  solutions  to  coupled  algebraic  Riccati-type  matrix  equations,  understood  as 
successive  approximations  to  the  original  equations  (69)-(70)  for  large  tf 

0  =  +  HXA  +  Q-  H1BR~1BtH1  +  nlH2BR~1BTH2 

+  n2[i3H2BR~1BTH3  +  [i2[i3H3BR~1BTH2  +  /d \H3BR~1BtH3 

0  =  AtH2  +  H2A  -  2^2H2BR-1BtH2  -  H2BR~1BtH1  -  H1BR~1BtH2 

-  fi3H2BR~1BTH3  -  n3H3BR~lBTH2  +  AH1GWGTH1 

0  =  AtH3  +  H3A  -  2[i3H3BR~1BTH3  -  H3BR~1BTH1  -  H1BR~1BTH3 

-  L i2H3BR~1BtH2  -  li2H2BR~1BtH3  +  6H1GWGtH2  +  6H2GWGTH1 . 

The  design  of  a  statistical  controller  using  first  three  cost  cumulants  was  carried  out  for  the  attitude  hold 
mode  when  /zi  =  1.0,  /i 2  =  2.0,  and  =  0.15.  Figure  2  illustrates  the  closed-loop  dynamical  behavior  of  the 
satellite  Euler  angles,  angular  rates  and  wheel  speeds  with  respect  to  time.  From  initial  values  as  indicated 
in  Figure  2,  the  satellite  transient  responses  went  to  zero  within  one  sixth  of  the  orbital  period.  Thus, 
the  statistical  control  method  performed  satisfactorily  at  the  attitude  hold  operation.  Figure  3  shows  the 
control  action  from  both  reaction  wheel  cluster  and  thrusters.  It  was  noted  that  the  action  of  the  statistical 
controller  was  large  at  the  start  and  settled  down  to  zero  very  fast. 

In  order  to  evaluate  the  attitude  hold  performance  of  the  satellite  against  space  debris  impact  and 
servicing  dock,  some  numerical  simulations  were  performed.  Figure  4  depicts  the  effects  of  two  space  debris 
collisions,  one  with  a  mass  of  0.5  kg  and  an  impact  speed  of  750  m/s  and  the  other  with  a  mass  of  0.3  kg 
and  an  impact  speed  of  800  m/s  on  the  Euler  angles,  angular  rates  and  wheel  speeds  of  the  satellite.  The 
alterations  in  the  satellite  pointing  angles  and  angular  rates  were  successfully  corrected  by  statistical  control 
activities  via  reaction  wheels  and  thrusters  as  can  be  seen  in  Figure  5.  And  the  statistical  controller  achieved 
finite-time  stabilization. 

Finally,  a  computer-aided  software  package  has  been  developed  to  provide  a  complete  statistical  descrip¬ 
tion  of  (43)  associated  with  the  Gauss-Markov  linear  dynamical  process  (41)  in  the  form  of  a  plot  of  the 
probability  density  function.  The  software  toolbox  helps  increase  the  depth  of  understanding  and  utility  of 
the  statistical  control  theory  in  terms  of  how  these  statistical  controllers  affect  the  overall  shape  of  cost  den¬ 
sities  after  control  selection  stages  are  taken  place.  As  illustrated  from  Figure  8,  the  shape  of  the  cost  density 
(43)  becomes  more  symmetric  as  the  skewness  weighting  /i3  increases  to  achieve  better  system  performance. 


VI.  Conclusions 

This  paper  has  addressed  the  application  question  by  applying  the  statistical  control  theory  to  the 
response  control  of  aerospace  structures  under  stochastic  and  impulsive  disturbances.  Although  the  statistical 
control  method  is  the  new  kid  on  the  block,  the  present  results  indicate  in  a  most  encouraging  way  that 
multiple  cost  cumulants  can  be  effectively  utilized  to  achieve  desirable  closed-loop  system  performance.  It 
is  hoped  that  these  early  results  serve  as  a  further  stimulus  for  this  very  interesting  approach. 
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Figure  2.  Closed-Loop  Responses  Due  To  Statistical  Controller  Design 
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Figure  3.  Reaction  Wheel  and  Thruster  Torques  Due  To  Statistical  Controller  Design 
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Figure  4.  Closed-Loop  Responses  Due  To  Space  Debris  of  0.5  kg  and  0.3  kg  and  Impact  Speeds  of  750  m/s 
and  800  m/s 
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Figure  5.  Control  Action  Due  To  Space  Debris  of  0.5  kg  and  0.3  kg  and  Impact  Speeds  of  750  m/s  and  800 
m/s 
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Figure  6.  Closed-Loop  Responses  Due  To  Servicing  Satellite  of  20  kg  and  Docking  Speed  of  5  m/s 
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Figure  7.  Control  Action  Due  To  Servicing  Satellite  of  20  kg  and  Docking  Speed  of  5  m/s 
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Figure  8.  Cost  Densities  Due  to  Statistical  Controllers 
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